使用GATK合并比较多个vcf文件

这不是简单的坐标合并，多个样本的vcf文件里面都只会记录对应样本的变异位点信息，但是因为样本检测手段不一致，可能是WGS或者WES，或者不会测序，意味着基因组某些区域可能是根本就没有被测序到，也就无从得知对应位点的碱基信息了。

这个时候的合并，要从bam文件开始，同一个坐标在某个样本既可能是野生型，杂合或者纯合突变，也有可能是该位点并没有任何reads覆盖。

如果不想自己写脚本去探究，那么gatk的HaplotypeCaller的GVCF模式正好能满足要求！

首先输出gvcf文件

代码如下，把所有的样本的bam文件都走一波gatk的HaplotypeCaller的GVCF模式，得到后缀为 _raw.g.vcf的文件（这个文件非常占硬盘空间哦）

1. GENOME=/home/jianmingzeng/reference/genome/human_g1k_v37/human_g1k_v37.fasta

2. TMPDIR=/home/jianmingzeng/tmp/software

3. GATK=/home/jianmingzeng/biosoft/GATK/GenomeAnalysisTK.jar

4. DBSNP=/home/jianmingzeng/annotation/variation/human/dbSNP/All_20160601.vcf.gz

5. ls /home/jianmingzeng/data/project/wes/bamFiles/*bam |while read id

6. do

7.    file=$(basename $id )

8.    sample=${file%%.*}

9.    echo $file $sample

10.    start=$(date +%s.%N)

11.    echo HaplotypeCaller `date`

12.    java -Djava.io.tmpdir=$TMPDIR   -Xmx40g -jar $GATK -T HaplotypeCaller  \

13.    -R $GENOME -I $id --dbsnp $DBSNP -ERC GVCF -gt_mode DISCOVERY \

14.    -stand_emit_conf 10 -o  ${sample}_raw.g.vcf

15.    echo HaplotypeCaller `date`

16.    dur=$(echo "$(date +%s.%N) - $start" | bc)

17.    printf "Execution time for HaplotypeCaller : %.6f seconds" $dur

18.    echo

19. done

理论上只要是bam文件里面表示该样本的基因组上面 某个位点被覆盖过，就会输出该位点的信息，无论其是否是突变。

这个输出的gvcf文件格式并不需要解释，也不需要理解，反正就是中间文件，当然，也欢迎有求知欲的同学继续深入了解哈。

合并多个样本的gvcf信息

因为每个样本的每个位置信息都被记录，所以这个时候的合并就很方便了，同样还是用gatk工具吧，代码如下：

1. GENOME=/home/jianmingzeng/reference/genome/human_g1k_v37/human_g1k_v37.fasta

2. TMPDIR=/home/jianmingzeng/tmp/software

3. GATK=/home/jianmingzeng/biosoft/GATK/GenomeAnalysisTK.jar

4. java -Djava.io.tmpdir=$TMPDIR   -Xmx40g -jar $GATK -T GenotypeGVCFs -nt 5  \

5. -R $GENOME   \

6. --variant  sample1.g.vcf   \

7. --variant  sample2.g.vcf   \

8. --variant  sample3.g.vcf   \

9. --variant  sample4.g.vcf \

10. -o output.vcf

另一种方法

以前我在直播我的基因组里面提到过，我的基因组是5条lane的独立fastq数据，期初我是先分开比对，然后把bam文件merge起来，结果发现自己在找变异的时候输出的vcf文件里面，每个lane都给出了基因型信息，也就是说根本就没有把这些lane当成是同一个样本。按照这个思路，我们可以把不同样本的bam文件合并，然后直接找变异，只有我们没有把样本信息抹掉，就仍然是独立样本独立基因型。

如果真正要合并不同的lane的数据为一个样本，需要用AddOrReplaceReadGroups来抹掉lane信息。

1. ls jmzeng_*.bam > files.bamlist

2. samtools merge -@ 5  -b files.bamlist  merged.bam

3. samtools index merged.bam

4. ### AddOrReplaceReadGroups ###

5. java -Djava.io.tmpdir=$TMPDIR    -Xmx40g -jar $PICARD AddOrReplaceReadGroups \

6. INPUT=${sample}.bam OUTPUT=${sample}_tmp.bam   RGID=jmzeng  RGLB=lib_all  RGPL=illumina RGPU=x10  RGSM=jmzeng

7. mv ${sample}_tmp.bam ${sample}.bam

8. samtools index ${sample}.bam

当然，vcf文件只是记录变异而已，变异一定要进行注释：

1. ls  /home/jianmingzeng/data/project/wes/variation/*_filtered_PASS.*.vcf |while read id

2. do

3. file=$(basename $id )

4. sample=${file%.*}

5. echo $sample

6. java -jar ~/biosoft/SnpEff/snpEff/snpEff.jar -i vcf GRCh37.75 $id >snpEFF_output/${sample}.snpEff.vcf

7. mv snpEff_summary.html ${sample}.snpEff_summary.html

8. ~/biosoft/ANNOVAR/annovar/convert2annovar.pl -format vcf4old  $id  >${sample}.annovar_input

9. ~/biosoft/ANNOVAR/annovar/annotate_variation.pl -buildver hg19 --geneanno --outfile annovar_output/${sample}_anno ${sample}.annovar_input  ~/biosoft/ANNOVAR/annovar/humandb/

10. done

更多gatk教程见：

本文只是合并了，至于如何比较，请听下回分解！